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We consider possible resolutions of singularities in a contracting anisotropic universe for a class of higher 
derivative gravity theories. We give evidence that for our models the big crunch singularity may be replaced by 
a nearly flat Minkowski-like phase before the Universe reenters an anisotropic expansion in a time-symmetric 
manner. 
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I. MOTIVATION AND RESULTS 

As a consequence of the singularity theorems of Penrose 
|Q]] and Hawking [2], it is clear that Einstein's theory of grav- 
ity is incomplete. To be more specific: a spacetime obeying 
Einstein's field equations and some reasonable conditions on 
causality and energy contains incomplete, inextendable time- 
like or null geodesies and is thus singular in the Schmidt 
sense. 

In general, the singularities predicted by these theorems are 
quite difficult to analyze, since the theorems do not predict 
their nature more concretely. But divergence of some curva- 
ture invariants seems to be a generic feature at any singular- 
ity in spacetime. In fact for singularities reached on timelike 
curves it can be proven (see yfl) that the Riemann tensor be- 
comes infinite. Thus effectively limiting the curvature seems 
to be a promising step on the way to a singularity free theory. 

On the other hand, one expects that a more fundamental the- 
ory of gravity avoids singular solutions completely. Based on 
these observations we are led to construct approximations to 
the low energy effective theory of full quantum gravity as min- 
imal deformations of Einstein's general relativity by higher 
derivative terms in such a way that at least the curvature sin- 
gularities are removed for generic solutions. A generic feature 
of low energy effective actions as limiting cases of fundamen- 
tal theories of gravity is the appearance of higher derivative 
terms as well as nonlocal terms. At present we will not con- 
sider the latter. Similar higher derivative corrections as we use 
here are motivated from the low energy effective action of the 
massless boson sector of string theory in |01 and fH]. 

In general relativity there are infinitely many curvature in- 
variants, e.g. 

which all have to be bounded in a nonsingular theory The 
idea of the limiting curvature hypothesis is now to impose an 
explicit bound on a finite subset of invariants only and make 
sure, that their limiting value uniquely determines a nonsin- 
gular spacetime. In this way we ensure that any spacetime, 
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that would be singular in Einstein gravity, is deformed into a 
nonsingular solution and henceforth all invariants are finite. 

In practice, the limiting of a suitable invariant / can be 
achieved through the introduction of a Lagrange multiplier 
field (f> with an appropriate potential V (</>). This idea was used 
by Brandenberger, Mukhanov, and Sornborger jg] for homo- 
geneous and isotropic spacetimes. 

In a Friedmann-Robertson- Walker cosmos with Hubble pa- 
rameter H = a I a, the invariant I 2 = AR^R^ - R 2 - 12H 2 
vanishes if and only if H — const, i.e. spacetime is de Sitter. 
More generally, I 2 — characterizes de Sitter spacetime as 
long as the Weyl tensor vanishes. Furthermore we remark that 
I2 is non-negative for a wide class of spacetimes (in particular, 
spherically symmetric ones). 

The detailed analysis [6] shows that the higher derivative 
theory with this invariant I2 has indeed no singular solutions, 
if we restrict ourselves to homogeneous and isotropic space- 
times. In this paper we investigate possible generalizations of 
this desingularization to anisotropic spacetimes. It is impossi- 
ble to describe the corresponding dynamical system by analyt- 
ical methods. For the model we consider in this work, we find 
evidence for the following generic behavior by a combination 
of analytical and numerical methods: (i) If the Lagrange mul- 
tiplier field <p grows large before the anisotropy becomes big 
then de Sitter spacetime is an asymptotic attractor in the far 
future, (ii) If the initial anisotropies are small then the lin- 
earized evolution equations are compatible with the attractor 
solution, (iii) If the initial anisotropies are of order one, then 
a generic solution passes through a phase which is approxi- 
mately Minkowski spacetime in finite time before returning 
to an anisotropic Universe in a time-symmetric fashion. This 
solution is nonanalytic. In particular, while all components of 
the Riemann tensor are finite, their derivatives are not. 

Apart from cosmological applications the higher derivative 
theory described here can also be applied to the interior of a 
black hole since the spacetime inside the horizon of a black 
hole resembles an anisotropic contracting, homogeneous uni- 
verse. The mechanism suggested here would suffice to build 
nonsingular black holes. The surrounding of the singularity 
would then be replaced by a region which approaches flat 
Minkowski spacetime interpolating between a black hole and 
a white hole solution. 
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II. THE MODEL 



An invariant which measures the anisotropy of a space- 
time is the squared Weyl tensor C 2 = C IJVTO -C f ' VTa ~ . To single 
out de Sitter spacetime from the class of homogeneous and 
anisotropic spacetimes, we can thus choose the sum of C 2 , 
whose vanishing guarantees isotropy, and the non-negative in- 
variant I2 selecting de Sitter from the isotropic solutions. We 
thus choose the higher derivative action to be of the form 



S = -j^ j(R-I<P + V(0)) y[=gd A x. 

where the invariant / is given by 

/ = h + 3C 2 = AR^R!" -R 2 + 3C MVTcr C l ™ rcr . 



(1) 



(2) 



Here I2 is the invariant chosen in the isotropic case. For an 
isotropic Universe the invariant / is, in fact, a perfect square. 
The factor 3 in front of the Weyl tensor squared is to simplify 
numerical factors, it does not alter the equations substantially, 
if one uses a different (but positive) factor. For the metric 
of an anisotropic and homogeneous cosmology (one special 
direction), we make the ansatz 

di 2 = -v(f) 2 df 2 + a(t) 2 dx 2 + /3{t) 2 [Ay 2 + dz 2 ) , (3) 

where the introduction of the lapse function v(t) is convenient 
to derive the field equations and set to 1 later. In this metric 
the invariant / takes the form 

/ = 4^(f + F 2 - H 2 f + 2 (H + H 2 - HF) 2 

+ {H-F-F 2 + HF) 2 ). (4) 

Here we introduced Hubble parameters H(t) - | and F(t) = ~ 
and set v(f) = 1 in the second line. Since / is the sum of 
three squares it vanishes if and only if all three terms vanish 
separately. There are two ways for this to happen. Either 



H(t) = F(t) = H = const.. 



(5a) 



which corresponds to the usual metric of spacially flat de Sit- 
ter spacetime (or Minkowski spacetime for H = 0), or 



H{t) = 0, 



Fit) = ? 



(5b) 



which describes Minkowski spacetime in unusual (singular) 
Kasner-like coordinates. 

Let us first discuss the conditions for approaching a de Sit- 
ter Universe. If Vi<p) rises slower than linear in for — > 00, 
the equations of motion derived from ([T]i will enforce de Sitter 
spacetime as soon as grows large. This is certainly achieved 
by a potential which approaches a constant value as goes 
to infinity. On the other hand, to recover the usual Einstein 
theory of gravity as a low curvature limit, after elimination of 
the unphysical field only higher order terms should remain. 
The auxiliary field is determined by the constraint equation 



Expanding V oc m + 0" 1+1 + ■ • • for small and keeping only 
the lowest order terms yields 



/ (R + AF n/< "'- l) ) yTgd 4 x for m > 1 
JiR + B + CI) V=gd 4 x for m = 1 



(7) 



with some constants A, B, C. Since / is itself quadratic in met- 
ric components, the correction term is of the desired order 
provided m > 1. The case m = 1 reduces to Einstein gravity 
with cosmological constant at low curvature. From the con- 
straint equation (O and the fact that the invariant / is strictly 
non-negative, it follows that V has to be a monotonically in- 
creasing function. 

A potential satisfying the conditions given above is 



Vi4>) = Hq 



4> sgn(0) 
(l + 



(8) 



More generally we may assume that the potential has an 
asymptotic power series in ^ of the form 



(9) 



V(d>) =m\\ + — r + 



where A is a positive constant and n > (the special case ([8]) 
corresponds to A = 1 and n = j). We shall see below that 
the existence of an asymptotic de Sitter attractor requires that 
n > 1 . On the other hand one expects that the appearance of 
fractional powers in the asymptotic expansion (O is a source 
for nonanalytic behavior. We will see below that this is indeed 
the case. 



III. FIELD EQUATIONS 

Variation of the higher derivative action ([1]) with respect to 
4> gives (after setting v = 1) 

(F + F 2 - H 2 f + 2 (H + H 2 - HF) 2 

2 1 

+ {F-H + F 2 -HF) = -V'i(p), (10a) 

while the 00-component of the Einstein equation takes the 
form 

= 8^ [iH + F)H + HF + 2HiH - F) 2 ) 

+40(2(7/ + F)H + 2HF - H 2 - 2HF + AHiH + 2F)H 
+ 2HiH + 2F)F -H 4 + 3H 2 F 2 - 2HF 3 ) 



-2H 1 - 4HF - V, 



(10b) 



/ = V'i<p). 



(6) 



where, again, we set v = 1 after variation. 

The equations arising from variation w.r.t. the spatial com- 
ponents of the metric a and B contain third derivatives of F 
and H but given ( llOal i and ( I lObb only one of the two spatial 
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equations is independent. Thus we can choose a linear combi- 
nation such that no third derivatives of F appear (after setting 
v(f) =1): 

=80 (-H + H(H - F)\ 

+40(-4# -2(2H + 3F) H - 2HF 

+ 4H 3 - 6FH 2 + 2F 2 H) 

+4<j>[-2H - 2 (5H + F)H- IB 2 - 2HF 

- 2 (8/Y 2 + 4HF -F 2 )ti-2 [h 2 - 3HF) F 



(10c) 



- 3H 4 - 2H 3 F + 5H 2 F 2 + 2HF 



') 



+4H + 6H £ + V. 



A closer look at the equations dTOb reveals that the total dif- 
ferential order is 5. It is possible to eliminate the Lagrange 
multiplier field with the help of ( II Oat , which is algebraic 
in 0. The resulting system of two differential equations for F 
and H (or F and H ) contains the inverse function of the po- 
tential. Since we do not want to specify the potential beyond 
its asymptotic properties we choose to keep in the equations. 

Let us now outline the strategy for solving the system of 
differential equations (ITOb . First we solve (II Oat for F. This 
involves taking the square root since ( II Oat is a quadratic func- 
tion in F. We then use the resulting equation to eliminate F in 
(llOat and solve for 0. Using these equations we can eliminate 
all time derivatives of F and in the spatial component of the 
field equation and thus obtain an equation for H. 



tuting this expansion into (llOab gives in leading order / = h. 1 
The second derivatives of the Hubble parameters can be ex- 
pressed through derivatives of 0: 



H 



" 2 



0(b+D/2 

J^r 



F = - 

With this expansion the 00-equation ( II Obi ) reduces to: 

1 J 1 \\ 



(15a) 
(15b) 



=4j>\H 1 ^-(j + 2h) 



0(«+l)/2 



+ 



+4M3H 2 {f + 2h) 



1 



0("+l)/2 
1 



+ 



03(»+l)/4 
1 



3, < 



n+l)/4 



(16) 



2 \0(" +1 >/ 4 



Neglecting higher order terms we get 



■ + H-(n + 3)0 - ^ + */V " +1)/2 = 0. 
2 ' Y f + 2h 



(17) 



For n < 1 a collapsing universe (// < 0) then leads to an 
asymptotically exponential growth 



• 3(ti + 3) , 
(p{t) oc exp — ■— — -Ht\, 



(18) 



IV. LARGE <f> LIMIT 

The case of large <p bears special interest since this limit 
corresponds to the regime where the invariant has to be small. 
We will discuss the two configurations ( 15at and d5bl with van- 
ishing invariant / separately. In the first case we will assume 
that at a given time spacetime is approximately de Sitter. As- 
suming the existence of a solution 0(f) we can replace the in- 
dependent variable t by 0, at least locally. The Hubble param- 
eters and their derivatives can now be expanded in powers of 
i. The 0-equation ( llOal i restricts the lowest order that can ap- 
pear in these expansions, since each square ist at most of the 
order of (4) 



Wl+l 



_ - h h 

H ~ H + 0(n+l)/4 + 0(«+l)/2 + ' 

f f 

F =H+ / ' + , „,„ + ■ 



,(«+l)/4 ^(«+l)/2 



F 



h 



0(n+l)/2 
/ 

0(n+l)/2 



(ID 
(12) 

(13) 

(14) 



These expansions are not independent, since H = tjf, which 
will be used later to determine the coefficients h and /. Substi- 



so that increases forever. For n — 1 there is a second term of 
order in ( fT7b . which leads to an asymptotic behavior of the 
form 



0(f) x exp 



3^ + ^o 2 3 g 

\mj+2h) 2 ) 



(19) 



This means for n = 1 = oo (i.e. de Sitter spacetime) is an 

attractor only if / + 2h < or / + 2h > \ + iS.. These 
conditions can however not be assumed a priori. They can 
only be verified once the asymptotic solution is known. This is 
in contrast to the isotopic model [6] where exponential growth 
of is generic. 

With the help of the asymptotic equation for 



-3/V0 + 0(0 3/4 ) 



(20) 



we can check the consistency of our ansatz ( fTTl i for H and H 
(and analogously for F and F), since 



H 



dH 



n -fh 



0(B+D/4 ^ 0( 



n+D/2 , 



+ 



1 



3 (' 



n+l)/4 



(21) 



Actually selfconsistency later requires h = 0. 
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Figure 1: Qualitative behavior of a possible solution bouncing 
through Minkowski 

spacetime. <p{i) diverges symmetrically around t = while F(t) 
diverges symmetrically. H{t) is smooth at t = 0. 



By comparing coefficients we get 
h = 



h=-(n+ l)Hh (22a) 



/=0 f=-(n+l)Hf (22b) 

therefore the leading terms in the expansion of H, F, H, F in 
(fTTT > are all determined in terms of h and /. 

Let us now turn to the second possibility (l5bl for vanishing 
invariant / at large 0. Here the Hubble parameters approach 
Minkowski spacetime in singular coordinates Oat , i.e. 



H( t) -> 0, 



F{tf + F(t) -> 0, 



fo- 



(23) 



Without restricting the generality we can assume fo = 0. In 
order to see if this behavior can be realized in a global solu- 
tion we are then looking for a solution of the field equations 
approaching H — and F — j for f — > 0. In this case the field 
equations immediately imply — > j. To continue we observe 
that the field equations are invariant under the discrete trans- 
formation simultaneously sending f — > -f, H(f) — > -H(t), 
F(t) — > ~F(t), and 0(f) — > 0(f). Approaching fo = from 
below with these asymptotics, we thus suggest the behavior 
sketched in Fig.Q] 

To make this discussion more precise let us rewrite the 
equations of motion in terms of G = l/F, H, and 0=1/0 
which are more suitable for analyzing the small f behavior. 
The proposed asymptotic solution is now G = t, = |f|, and 
H = 0. Therefore we seek approximate solutions of the form 



G - t + ai t 1+ai + a 2 t l+ai+ai + a 3 t l+ai+a ^ +..., 
H = b ] t j3 < + bof'^ 2 + b^' +/3 ^ + ..., 



= +t + C\t 



i+n 



:C 2 t 



i+n +72 



:C3f i+y, +n+y, + 



(24) 
(25) 
(26) 



where a,, /J,-, and y, are positive constants and the two signs 
correspond to f % 0. In the following we will only discuss the 
case t > 0. For convenience we choose a specific potential 
with asymptotic expansion 

4V 2 10V 3 

01/2 + ~ 0372 



V(0) = Vj 



20 y 4 35 y 5 

+ IT- + 



(27) 



The potential V(0) = corresponds to the choice 

V, = Hq. Plugging this ansatz into the equations of motion 
and neglecting all obviously subleading terms we remain with 



W 



= — ? 3/2 + ft? (3 - 2/3 1 + 3fi)t 2(/Sl 

+ 2a l b x (\ +a x ){\ + Pi)r 
= V 1 + Sb 1 (2+J3 l -f 1 )/'- 4 , 
= y 1 -8f> 1 (6+/3 1 -4^+ jS 3 )^ 



+ 2a](\+a l yt 2 ' 



(28a) 

,2a, -4 

(28b) 
(28c) 



In ( I28bb and ( I28cb we can either choose ySi =4 and cancel 
both terms against each other or we have to choose fi\ < 4 as 
a zero of the two coefficients. We are left with two admissible 
solutions obeying f}\ > 0: 



a\ = — , Pi= 2, fli = ±^". 
11 

a\ P\= 4, fli = ±^3-. 



or 



(29a) 



bi = -. (29b) 



Proceeding in this way we obtain a perturbative solution such 
that the equations of motion are satisfied up to linear order at 
f = with 



G = t-^f'\ 
llbi 



H -b\t + (48/>iCi + Vi)r + 



= ±t + c\ f ± C2/ 1 for f % 0. 



4 2(1155/>ic 2 -38y 2 ) 



4235 



(30a) 

.9/2 

(30b) 
(30c) 



We have not been able to extend the solution as a power series 
in fractional powers of f beyond this order suggesting that, if a 
solution exists for a finite range of f, nonanalytic behavior of a 
different kind will be required. Note also that all components 
of the Riemann tensor (but not all its derivatives) are finite at 
f = 0. 

To summarize, we found that if is already large while the 
two Hubble parameters are comparable the proposed limiting 
curvature procedure works perfectly well. In the case where 
the leading order correction to the potential has an exponent 
n < 1, large implies that spacetime is nearly de Sitter and 
the evolution of extends to infinite future while growing ex- 
ponentially and thus forcing spacetime to approach de Sitter 
even more. Because of the exponential growth of it takes an 
infinite time to actually reach the de Sitter end stage. On the 
other hand, if the difference of the two Hubble parameters is 
of order one while is large, the metric may approach Min- 
kowski spacetime in singular coordinates in finite time. Al- 
though this solution, if it exists, is nonanalytic at f = 0, it can 
be continued symmetrically through the "Minkowski phase" 
at f = 0. In Sec. [VI] we will further analyze the existence of a 
global solution using numerical methods. 



V. PHASE SPACE ANALYSIS 

In this section we discuss the global properties of the so- 
lutions of the equation of motion by means of a phase space 
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analysis. The equations of motion ( fTOb form a system of or- 
dinary differential equations with total differential order of 5. 
We eliminate derivatives of F from the equations (II Obi ) and 
( II Pel ) with the help of dlOab and all derivatives of (f> from 
Eq. JlOct by ( II Obi ). As usual we rewrite the equations into 
a system of first order differential equations for the vector 
u := F, H, H, H) and obtain 



it = V(u) 



(31) 



Since ( llOab contains F 2 solving the equations involves a 
square root of 



B 



V'(<p)-10F 2 H 2 +20FH 3 - 



\QH 4 + \2H{F-H)H- 10H 2 . 

(32) 

The reality of the solution is thus not guaranteed a priori. 
However, it turns out that the hypersurface B — is left in- 
variant under the flow of the vectorfield V, since 



du l 



= 0. 



(33) 



Hence, if we restrict the initial conditions it' to the domain 
of B > 0, i.e. real vector field V(W(0)), the solution will remain 
real for all times. 

Let us now discuss the existence of singular points. Crit- 
ical points where V(u) vanishes are problematic if they cor- 
respond to fixed points away from the asymptotic isotropic 
regime. Since V 3 = u 4 and V 4 = u 5 , we can reduce the prob- 
lem of searching for critical points V — to two dimensions 
by replacing u 4 and u 5 with 0, and u 2 by one of the four so- 
lutions of the equation V 2 = 0. It can be shown that the re- 
maining two functions V l (u l ,u 2 ) and V 5 (u l , u 2 ) do not vanish 
simultaneously. But V 5 (u l ,u 2 ) has a one-dimensional set of 
fractional zeros and a one- dimensional set of fractional poles. 
The poles typically correspond to branching points of the solu- 
tion. If more than one real solution meet at a branching point, 
then the Cauchy problem is not well defined. If furthermore 
the set of poles and zeros intersect then further complications 
arise since the details will now depend on how the solution 
approaches the singularity. We can not exclude the existence 
of such points; but a numerical investigation of | V 5 (m', m 2 )| in- 
dicates that whether the intersection is empty or not depends 
on the details of the potential V{<p). 



VI. NUMERICAL SOLUTIONS 

In this section we supplement the analytic discussion of the 
set of solutions to the higher derivative action by a numerical 
investigation. Of particular interest is the question of whether 
the nonanalytic local solution passing through the Minkow- 
ski stage can be embedded in a global solution and to deter- 
mine the fate of a generic anisotropic Universe. While we 
are not able to give a conclusive answer to these questions, 
we find that the numerical analysis supports the idea that a 
generic anisotropic initial condition will go through a Min- 
kowski phase before returning to the anisotropic Universe in 
a time-symmetric fashion. 



H = F 



Figure 2: Two numerical solutions (solid and dashed lines; H = F 
(thick, negative axis), cf> (thin, positive axis)) starting with isotropic 
initial conditions. 




Figure 3: Two numerical solutions (solid and dashed lines; H, F 
(thick), and H - F (thin) where the latter is rescaled with a factor of 
10) starting with initial conditions slightly perturbed from isotropy. 



Let us begin however with isotropic initial conditions, the 
behavior of the solutions of |0] is recovered — the Hubble 
parameter H{t) = F(t) approaches a constant while (f> grows 
exponentially. The two numerical solutions plotted in Fig. [2] 
were obtained using initial conditions 



HQS) = F(0) 



-0.1 
-0.3 



0(0) = 0.1. 



For these and further numerical solutions we chose the poten- 
tial to be of the form ([8]l with the constant set to Ho = 10. 

Initial conditions prescribing a slightly perturbed isotropic 
universe evolve to an isotropic solution — the perturbations 
are damped and the Hubble parameters and <p show the same 
asymptotic behavior as in the isotropic case. Fig.[3]shows the 
difference H - F decay while the two Hubble parameters ap- 
proach a constant. The unphysical field <p is not plotted since 
its behavior is similar to the isotropic case. The initial condi- 
tions used were 



ff(0) = 



1-0.1 
1-0.3 



FQS) = 



1-0.12 
I -0.28 



0(0) = 0.1. 
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This numerical prediction can be supported by consider- 
ing the linearization of the system of equations around an 
isotropic background with F(t) = H(t) satisfying the 00- 
component of the field equations with <p given by the con- 
straint equation, we obtain a system of quasilinear equations. 
The linearized system of field equations decouples after intro- 
duction of the new dependent variables x(i) = 6F(t) — 6H(t) 
and y(t) = 6F(t) + 26H(t). The equation for the linearized 
anisotropy x(f) is 

= jc + 2\3H + - )x + \15H 2 + 6H + 3H- + - - — )x 

[ <t>) \ <?> <f> 2<p) 

+ 3 (6H 3 + 1HH + ti) x. (34) 

Hence, a hint on the linear stability of the isotropic solution 
can be derived from the eigenvalues of this differential equa- 
tion. With the help of the background field equations, we can 
express the eigenvalues in terms of H and <j>. In a contract- 
ing universe (H < 0) the real part of at least one eigenvalue is 
positive if and/or H are small which might result in an insta- 
bility. Hence, we conclude that the isotropic solution is stable 
if <p is large. Furthermore, a numerical analysis shows that 
the linear perturbations stay small long enough for the back- 
ground field 4> to become large. Note that, since the square 
of the Weyl tensor vanishes in the isotropic background to- 
gether with its first derivative, there are no contributions from 
the Weyl tensor in the linearized equations. Thus, they are 
no longer applicable, as soon as the background field (f> be- 
comes large and thus the nonlinear constraint equation ( II Oat 
becomes important. 

The numerical evolution of generic initial conditions breaks 
down at a pole like singularity in F, while H approaches 0, 
and (p has just passed a minimum and thus seems to be ap- 
proximately constant. But from (1 1 Oat we imply that <p oc ± 
as soon as F approaches the singularity further. This behav- 
ior resembles the aforementioned singular coordinate system 
of Minkowski spacetime. We therefore suggest that these nu- 
merical solutions could be patched to the approximate solu- 
tion (f30b describing a kind of bounce through flat Minkowski 
spacetime. 

In Fig. 2] this continuation is sketched. The numerical so- 
lution is integrated from the left until it stops. Using the free 
parameters that are left in ( f30l > the asymptotical solution is fit- 
ted approximately. We do not expect to be able to patch the 
two sets of functions together smoothly since there is no over- 
lap between the range of validity of the numerical solution 
and the approximate solution ( f30b at t — to- In particular H 
obviously does not meet. Nevertheless the picture hints on the 
existence of a continued solution. 



VII. CONCLUSIONS 

In this paper we considered a higher derivative theory of 
gravity in which the square of the Weyl tensor together with 
another invariant, h, whose vanishing guarantees de Sitter 
spacetime, is bounded with the help of a Lagrange multiplier 



H 

F 




t 



Figure 4: Sketch of the proposed matching of the numerical solution 
on the left and the possible approximate continuation on the right. 
There is no overlap between the respective regions of validity of the 
two approximations. 

field. For a generic model of this type we find two qualita- 
tively different end points of a contracting, anisotropic Uni- 
verse: If the initial anisotropies are small then the anisotropies 
stay small in the linearized regime for long enough allowing 
the Lagrange multiplier field to grow. We found convincing 
evidence that for this set of initial configurations the non- 
linear effects will then suppress the anisotropies leading to 
an asymptotic contracting de Sitter spacetime for late times. 
Since the Weyl tensor appears only quadratically in the higher 
derivative theory the anisotropies are not suppressed in the 
linear regime. A possible way to improve on this in the linear 
regime might be to introduce separate Lagrange multipliers 
for the Weyl tensor and I2, respectively. This would probably 
help to enlarge the regions of applicability of the linearized 
theory. 

On the other hand, if the initial anisotropies are of order 
1, then the linearized theory does not apply. A combination 
of analytical methods and numerical approximations in this 
regime produced circumstantial evidence that a global solu- 
tion may exist in this case which interpolates between a con- 
tracting anisotropic Universe and a nearly Minkowski phase 
in finite time before returning to an anisotropic spacetime in 
a time-symmetric fashion. While we are not able to prove 
the existence of such a solution globally at present, we think 
that a better understanding of this solution, if it exists, could 
have important applications in the resolution of cosmologi- 
cal singularities as an alternative to the so-called bounce so- 
lutions. Furthermore the same mechanism could be used to 
resolve the spacelike singularity inside a Schwarzschild black 
hole. In our model we find evidence that the surrounding of 
the singularity is replaced by a region which approaches flat 
Minkowski spacetime interpolating between a black hole and 
a white hole solution. This resembles the bridge between two 
large classical solutions obtained from a quantum geometri- 
cal discussion in 10]. In an alternative scenario proposed in 
|8y, the black hole singularity is replaced by a region which 
is approximately de Sitter spacetime. While we do not find 
evidence for this scenario in our model we can not exclude it 
at present. 
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